## The code below replicates all tables and figures 
## (except Figure 1; see Replication_AB.R) 
## from the main manuscript and appendix in
## Yuree Noh and Marwa Shalaby, "Who Supports Gender Quotas in Transitioning and Authoritarian States in the Middle East and North Africa?" Comparative Political Studies

############################
## Load Data and Packages ##
############################

rm(list=ls())
library(foreign)
mor = read.csv("morocco_clean_publish.csv")
tun = read.csv("tunisia_clean_publish.csv")

library(stargazer)
library(ggplot2)
library(MASS)
library(dplyr)
library(tidyr)
library(mlogit)
library(dotwhisker)
library(broom)
library(reshape)
library(nnet)


# Table 1
t1 = table(tun$quota)
t1a = cbind(t1, prop.table(t1))
t2 = table(mor$quota)
t2a = cbind(t2, prop.table(t2))
library(xtable)
xtable(cbind(t1a, t2a))




# Figure 2
p1 = cbind(c(0,0,0,0,1,1,1,1), c("No quota","Less","As is", "More", "No quota","Less","As is", "More"), c(11.72,18.38,40.40,29.49,5.47,5.67,28.95,59.92))
colnames(p1) = c("Gender", "Support", "value")
p1 = as.data.frame(p1)
p1[,3] <- as.numeric(as.character(p1[,3]))

p1$Support = factor(p1$Support, levels=c("No quota","Less","As is", "More"))
pp1 = ggplot(p1, aes(x=Support, y=value, fill=Gender)) + geom_bar(stat="identity", position="dodge") + scale_fill_manual(name="Gender", labels=c("Male", "Female"), values =c("grey66", "grey20")) + xlab("") + ylab("%") + geom_text(size=2.8, aes(label=paste0(sprintf("%0.f", value), "%")), position = position_dodge(width = 0.9), vjust=-0.3, size=3.5) + ggtitle("Morocco") + scale_y_continuous(limits = c(0,62)) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) 
pp1

p2 = cbind(c(0,0,0,0,1,1,1,1), c("No quota","Less","As is", "More", "No quota","Less","As is", "More"), c(25.37,28.36,41.56,4.71,25.34,17.59,45.76,11.32))
colnames(p2) = c("Gender", "Support", "value")
p2 = as.data.frame(p2)
p2[,3] <- as.numeric(as.character(p2[,3]))

p2$Support = factor(p2$Support, levels=c("No quota","Less","As is", "More"))
pp2 = ggplot(p2, aes(x=Support, y=value, fill=Gender)) + geom_bar(stat="identity", position="dodge") + scale_fill_manual(name="Gender", labels=c("Male", "Female"), values =c("grey66", "grey20")) + xlab("") + ylab("") + geom_text(size=2.8, aes(label=paste0(round(value,0), "%")), position = position_dodge(width = 0.9), vjust=-0.3, size=3.5) + ggtitle("Tunisia") + scale_y_continuous(limits = c(0,62)) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) 
pp2

library(ggpubr)
ggarrange(pp1, pp2, common.legend = TRUE, legend = "right") 




# Table 2
mor$quota_dummy = ifelse(mor$quota > 2, 1, 0)
tun$quota_dummy = ifelse(tun$quota > 2, 1, 0)
table(mor$quota_dummy)
table(tun$quota_dummy)
mod1 = glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age , family = binomial(), data=tun)
mod2 = glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda , family = binomial(), data=tun)
mod3 <- glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age , family = binomial(), data=mor)
mod4 <- glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam , family = binomial(), data=mor)
# log odds
# source: https://github.com/cimentadaj/cimentadaj/blob/master/R/stargazer2.R
stargazer2 <- function(model, odd.ratio = FALSE, ...) {
  if(!("list" %in% class(model))) model <- list(model)
  
  if (odd.ratio) {
    coefOR2 <- lapply(model, function(x) exp(stats::coefficients(x)))
    seOR2 <- lapply(model, function(x) exp(stats::coefficients(x)) * summary(x)$coef[, 2])
    p2 <- lapply(model, function(x) summary(x)$coefficients[, 4])
    stargazer::stargazer(model, coef = coefOR2, se = seOR2, p = p2, ...)
    
  } else {
    stargazer::stargazer(model, ...)
  }
}
stargazer2(mod1, odd.ratio=TRUE)
stargazer2(mod2, odd.ratio=TRUE)
stargazer2(mod3, odd.ratio=TRUE)
stargazer2(mod4, odd.ratio=TRUE)
# combined by hand





# Figure 3 & Figure 4
library(sjPlot)
fig1 = get_model_data(mod2, type = "pred", terms = c("gov_index"))
plot.fig1 = ggplot(fig1, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) + 
  geom_line(aes(x = x, y = predicted, group = group), 
            lwd = 1/2, colour="black") + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              alpha=0.1, colour = "black", linetype = 2) +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Tunisia") + 
  xlab("Government Performance") + ylab("Predicted Probabilities")

fig2 = get_model_data(mod4, type = "pred", terms = c("gov_index"))
plot.fig2 = ggplot(fig2, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) + 
  geom_line(aes(x = x, y = predicted, group = group), 
            lwd = 1/2, colour="black") + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              alpha=0.1, colour = "black", linetype = 2) +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Morocco") + 
  xlab("Government Performance") + ylab("Predicted Probabilities")

fig3 = get_model_data(mod2, type = "pred", terms = c("elec_free"))
plot.fig3 = ggplot(fig3, aes(ymin = 0.3, ymax = 1)) + 
  geom_line(aes(x = x, y = predicted, group = group), 
            lwd = 1/2, colour="black") + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              alpha=0.1, colour = "black", linetype = 2) +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Tunisia") + 
  xlab("Trust in Elections") + ylab("Predicted Probabilities")

fig4 = get_model_data(mod4, type = "pred", terms = c("elec_free"))
plot.fig4 = ggplot(fig4, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) + 
  geom_line(aes(x = x, y = predicted, group = group), 
            lwd = 1/2, colour="black") + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              alpha=0.1, colour = "black", linetype = 2) +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Morocco") + 
  xlab("Trust in Elections") + ylab("Predicted Probabilities")

library(ggpubr)
ggarrange(plot.fig1, plot.fig2,
          ncol = 2)
ggarrange(plot.fig3, plot.fig4,
          ncol = 2)






# Table 3

mod5 = glm(quota_dummy ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda , family = binomial(), data=tun, subset = (female == 1))
summary(mod5)
mod6 <- glm(quota_dummy ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda , family = binomial(), data=tun, subset = (female == 0))
summary(mod6)
mod7 = glm(quota_dummy ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam , family = binomial(), data=mor, subset = (female == 1))
summary(mod7)
mod8 <- glm(quota_dummy ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam , family = binomial(), data=mor, subset = (female == 0))
summary(mod8)

stargazer2(mod5, odd.ratio=TRUE)
stargazer2(mod6, odd.ratio=TRUE)
stargazer2(mod7, odd.ratio=TRUE)
stargazer2(mod8, odd.ratio=TRUE)
# combined by hand







# Figure 5 & Figure 6
fig5 = get_model_data(mod5, type = "pred", terms = c("gov_index")) # tun gov women 
fig5 = fig5[c(1,4,7,10),]
fig6 = get_model_data(mod6, type = "pred", terms = c("gov_index")) # tun gov men
fig6 = fig6[c(1,4,7,10),]
fig7 = get_model_data(mod7, type = "pred", terms = c("gov_index")) # mor gov women 
fig7 = fig7[c(1,4,7,10),]
fig8 = get_model_data(mod8, type = "pred", terms = c("gov_index")) # mor gov men
fig8 = fig8[c(1,4,7,10),]
fig5a = get_model_data(mod5, type = "pred", terms = c("elec_free")) # tun elec women 
fig6a = get_model_data(mod6, type = "pred", terms = c("elec_free")) # tun elec men 
fig7a = get_model_data(mod7, type = "pred", terms = c("elec_free")) # mor elec women 
fig8a = get_model_data(mod8, type = "pred", terms = c("elec_free")) # mor elec women 

# tunisia gov
plot1 = ggplot(fig5, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) + 
  geom_line(aes(x = fig5$x, y = fig5$predicted), 
            lwd = 1/2, colour="black", linetype="longdash")  + geom_point(aes(x = fig5$x, y = fig5$predicted), shape=87, size=4, color = "black") +
  geom_line(aes(x = fig6$x, y = fig6$predicted), 
            lwd = 1/2, colour="black", linetype="dotted") +  geom_point(aes(x = fig6$x, y = fig6$predicted), shape=77, size=4, color = "black") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Tunisia") + 
  xlab("Government Performance") + ylab("Predicted Probabilities")

# morocco gov
plot2 = ggplot(fig7, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) + 
  geom_line(aes(x = fig7$x, y = fig7$predicted), 
            lwd = 1/2, colour="black", linetype="longdash") + geom_point(aes(x = fig7$x, y = fig7$predicted), shape=87, size=4, color = "black") +
  geom_line(aes(x = fig8$x, y = fig8$predicted), 
            lwd = 1/2, colour="black", linetype="dotted") + 
  geom_point(aes(x = fig8$x, y = fig8$predicted), shape=77, size=4, color = "black") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Morocco") + 
  xlab("Government Performance") + ylab("Predicted Probabilities")

# tunisia elec
plot3 = ggplot(fig5a, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) +
  geom_line(aes(x = fig5a$x, y = fig5a$predicted), 
            lwd = 1/2, colour="black", linetype="longdash")  + geom_point(aes(x = fig5a$x, y = fig5a$predicted), shape=87, size=4, color = "black") +
  geom_line(aes(x = fig6a$x, y = fig6a$predicted), 
            lwd = 1/2, colour="black", linetype="dotted") +  geom_point(aes(x = fig6a$x, y = fig6a$predicted), shape=77, size=4, color = "black") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Tunisia") + 
  xlab("Trust in Elections") + ylab("Predicted Probabilities")

# morocco elec
plot4 = ggplot(fig7a, aes(ymin = 0.3, ymax = 1, colour = as.factor(group))) + 
  geom_line(aes(x = fig7a$x, y = fig7a$predicted), 
            lwd = 1/2, colour="black", linetype="longdash")  + geom_point(aes(x = fig7a$x, y = fig7a$predicted), shape=87, size=4, color = "black") +
  geom_line(aes(x = fig8a$x, y = fig8a$predicted), 
            lwd = 1/2, colour="black", linetype="dotted") +  geom_point(aes(x = fig8a$x, y = fig8a$predicted), shape=77, size=4, color = "black") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "bottom") + 
  theme(legend.position="") + ggtitle("Morocco") + 
  xlab("Trust in Elections") + ylab("Predicted Probabilities")

ggarrange(plot1, plot2,
          ncol = 2)
ggarrange(plot3, plot4,
          ncol = 2)







# Figure A.1 and A.2

par(mfrow=c(1,4))
hist(mor$age, main="", xlab="Age", breaks=8, col="darkgrey")
abline(v = mean(mor$age), lty = 2)

table(mor$edu)
table1 <- matrix(c(21.29, 19.87, 33.03, 25.81), ncol=4,byrow=TRUE)
colnames(table1) <- c("None", "Primary", "Secondary", "Tertiary")
rownames(table1) <- ""
barplot(table1, ylim=c(0,100), ylab="Percentage", xlab="Highest Level of Education", col="darkgrey")

table(mor$married)
table1 <- matrix(c(51.75, 48.25),ncol=2,byrow=TRUE)
colnames(table1) <- c("Married","Other")
rownames(table1) <- ""
barplot(table1, ylim=c(0,100), ylab="Percentage", xlab="Marital Status", col="darkgrey")

table(mor$urban)
table1 <- matrix(c(58.67, 41.33),ncol=2,byrow=TRUE)
colnames(table1) <- c("Urban","Rural")
rownames(table1) <- ""
barplot(table1, ylim=c(0,100), ylab="Percentage", xlab="Urban/Rural", col="darkgrey")




par(mfrow=c(1,4))
hist(tun$age, main="", xlab="Age", breaks=8, col="darkgrey")
abline(v = mean(tun$age), lty = 2)

table(tun$edu)
table1 <- matrix(c(14.11, 26.81, 40.54, 18.53), ncol=4,byrow=TRUE)
colnames(table1) <- c("None", "Primary", "Secondary", "Tertiary")
rownames(table1) <- ""
barplot(table1, ylim=c(0,100), ylab="Percentage", xlab="Highest Level of Education", col="darkgrey")

table(tun$married)
table1 <- matrix(c(40.60, 59.40),ncol=2,byrow=TRUE)
colnames(table1) <- c("Married","Other")
rownames(table1) <- ""
barplot(table1, ylim=c(0,100), ylab="Percentage", xlab="Marital Status", col="darkgrey")

table(tun$urban)
table1 <- matrix(c(31.98, 68.02),ncol=2,byrow=TRUE)
colnames(table1) <- c("Urban","Rural")
rownames(table1) <- ""
barplot(table1, ylim=c(0,100), ylab="Percentage", xlab="Urban/Rural", col="darkgrey")



# Figures A.3-A.6
par(mfrow=c(1,1))

## Figure A.3
tuntablef = prop.table(table(tun$quota[tun$female == 1]))*100
tuntablef
# 25.33825 17.58918 45.75646 11.31611 
tuntablef2 = barplot(tuntablef, main = "", names.arg=c("No quota","Less than 50%","Just 50%", "More than 50%"), ylim=c(0,70), ylab="%")
text(x=tuntablef2, y=tuntablef, label=c("25.34%", "17.59%", "45.76%", "11.32%"), pos=3, cex=0.8, col="blue")

## Figure A.4
tuntablem = prop.table(table(tun$quota[tun$female == 0]))*100
tuntablem
# 25.373134 28.358209 41.561424  4.707233
tuntablem2 = barplot(tuntablem, main = "", names.arg=c("No quota","Less than 50%","Just 50%", "More than 50%"), ylim=c(0,70), ylab="%")
text(x=tuntablem2, y=tuntablem, label=c("25.37%", "28.36%", "41.56%", "4.71%"), pos=3, cex=0.8, col="blue")

## Figure A.5
mortablef = prop.table(table(mor$quota[mor$female == 1]))*100
mortablef
# 5.465587  5.668016 28.947368 59.919028
mortablef2 = barplot(mortablef, main = "", names.arg=c("No quota","Less than 60","Just 60", "More than 60"), ylim=c(0,70), ylab="%")
text(x=mortablef2, y=mortablef, label=c("5.47%", "5.67%", "28.95%", "59.92%"), pos=3, cex=0.8, col="blue")

## Figure A.6
mortablem = prop.table(table(mor$quota[mor$female == 0]))*100
mortablem
# 11.71717 18.38384 40.40404 29.49495 
mortablem2 = barplot(mortablem, main = "", names.arg=c("No quota","Less than 60","Just 60", "More than 60"), ylim=c(0,70), ylab="%")
text(x=mortablem2, y=mortablem, label=c("11.72%", "18.38%", "40.40%", "29.49%"), pos=3, cex=0.8, col="blue")


## Table A.1
main = data.frame(
  
  GovernmentPerformance = mor$gov_index,
  TrustinElection = mor$elec_free,
  Female = mor$female,
  GenderEgalitarianism = mor$gender_index_soc,
  Islamist = mor$islamist, 
  Education = mor$edu,
  Urban = mor$urban,
  HouseholdFinance = mor$household_fin,
  Married = mor$married,
  Age = mor$age,
  PJD = mor$pjd,
  PAM = mor$pam,
  FemaleIntrvr = mor$intrvwr_f)

library(stargazer); stargazer(main, title="Descriptive Statistics (Morocco)", omit.summary.stat = c("p25","p75"), digits=1)


## Table A.2
main = data.frame(
  GovernmentPerformance = tun$gov_index,
  TrustinElection = tun$elec_free,
  Female = tun$female,
  GenderEgalitarianism = tun$gender_index_soc,
  Islamist = tun$islamist, 
  Education = tun$edu,
  Urban = tun$urban,
  HouseholdFinance = tun$household_fin,
  Married = tun$married,
  Age = tun$age,
  NidaaTounes = tun$nidaa,
  Ennahda = tun$ennahda,
  FemaleIntrvr = tun$intrvwr_f)

library(stargazer); stargazer(main, title="Descriptive Statistics (Tunisia)", omit.summary.stat = c("p25","p75"), digits=1)




### Tables A.3 - A.6
cor.test(tun$gov_index, tun$female, use = "complete.obs")
fit1 <- lm(gov_index ~ female, data = tun)
summary(fit1)

cor.test(tun$elec_free, tun$female, use = "complete.obs")
fit2 <- lm(elec_free ~ female, data = tun)
summary(fit2)

cor.test(tun$gov_index, tun$elec_free, use = "complete.obs")
fit3 <- lm(gov_index ~ elec_free, data = tun)
summary(fit3)

cor.test(mor$gov_index, mor$female, use = "complete.obs")
fit4 <- lm(gov_index ~ female, data = mor)
summary(fit4)

cor.test(mor$elec_free, mor$female, use = "complete.obs")
fit5 <- lm(elec_free ~ female, data = mor)
summary(fit5)

cor.test(mor$gov_index, mor$elec_free, use = "complete.obs")
fit6 <- lm(gov_index ~ elec_free, data = mor)
summary(fit6)

# Table A.7
stargazer(fit1,fit2,fit3)
# Table A.8
stargazer(fit4,fit5,fit6)





## Table A.9
stargazer(mod1, mod2, mod3, mod4,
          table.placement="H", omit.stat=c("f", "ser", "rsq", "adj.rsq"),
          dep.var.labels.include = F,
          covariate.labels = c("Government \nPerformance", "Trust in \nElections", "Female", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Age", "Nidaa \nTounes", "Ennahda", "PJD", "PAM"), no.space = F)


## Table A.10
stargazer(mod5, mod6, mod7, mod8,  
          table.placement="H", omit.stat=c("f", "ser", "rsq", "adj.rsq"),
          dep.var.labels.include = F,
          covariate.labels = c("Government \nPerformance", "Trust in \nElections", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Age", "Nidaa \nTounes", "Ennahda", "PJD", "PAM"), no.space = F)

## Table A.11
# pooled results - merge data
mor$country = "mor" # morocco dummy
tun$country = "tun"
library(dplyr)
dat = bind_rows(mor,tun)
mod_pool = glm(quota_dummy ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + female + age + intrvwr_f + factor(country), family = binomial(), data=dat)
stargazer(mod_pool, table.placement="H", omit.stat=c("f", "ser", "rsq", "adj.rsq"),
          dep.var.labels.include = F, 
          covariate.labels = c("Government \nPerformance", "Trust in \nElections", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Female", "Age", "Female \nInterviewer", "Tunisia"), no.space = F)



## Table A.12
### robustness - govt index with indiv questions

mod1a = glm(quota_dummy ~ gov_handle + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda, family = binomial(), data=tun)
mod1b = glm(quota_dummy ~ gov_corruption + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda, family = binomial(), data=tun)
mod1c = glm(quota_dummy ~ gov_response + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda, family = binomial(), data=tun)
mod2a = glm(quota_dummy ~ gov_handle + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam, family = binomial(), data=mor)
mod2b = glm(quota_dummy ~ gov_corruption + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam, family = binomial(), data=mor)
mod2c = glm(quota_dummy ~ gov_response + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam, family = binomial(), data=mor)

stargazer(mod1a, mod1b, mod1c, mod2a, mod2b, mod2c, 
          table.placement="H", omit.stat=c("f", "ser", "rsq", "adj.rsq"),
          dep.var.labels.include = F,
          covariate.labels = c("Government \nHandle", "Government \nCorruption", "Government \nResponsive", "Trust in \nElections", "Female", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Age", "Nidaa \nTounes", "Ennahda", "PJD", "PAM"), no.space = F)





# Table A.13
#multinomial logit analysis with 3 categories
#### robustness
library(nnet)
mor$quota1 <- factor(mor$quota, ordered = FALSE )
mor$quota1 = relevel(mor$quota1, ref = "3")
tun$quota1 <- factor(tun$quota, ordered = FALSE )
tun$quota1 = relevel(tun$quota1, ref = "3")

mn1 = multinom(quota1 ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + female + age + intrvwr_f , data=tun)
summary(mn1)

mn2 = multinom(quota1 ~ gov_index + elec_free + gender_index_soc + islamist + edu + urban + household_fin + married + female + age + intrvwr_f, data=mor)
summary(mn2)

mn1.rrr = exp(coef(mn1))
mn2.rrr = exp(coef(mn2))
exp(coef(mn1))
exp(coef(mn2))

n1 = nrow(mn1$residuals)
n2 = nrow(mn2$residuals)
stargazer(mn1, mn2, coef=list(mn1.rrr, mn2.rrr), p.auto=FALSE,
          type="latex", 
          add.lines = list(c("\\textit{$N$}", n1,n1,n1,n2,n2,n2)),
          covariate.labels = c("Government \nPerformance", "Trust in \nElections", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Female", "Age", "Female \nInterviewer"), no.space = F)




# Table A.14
table(mor$A8_1)

# Table A.15
table(tun$prob1)

# Table A.16
# With Issue Controls
mod_issue_t = glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda + issue_f, family = binomial(), data=tun)
mod_issue_m <- glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam + issue_f, family = binomial(), data=mor)
stargazer(mod_issue_t, mod_issue_m,
          table.placement="H", omit.stat=c("f", "ser", "rsq", "adj.rsq"),
          dep.var.labels.include = F,
          covariate.labels = c("Government \nPerformance", "Trust in \nElections", "Female", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Age", "Nidaa \nTounes", "Ennahda", "PJD", "PAM", "Feminine \nIssues"), no.space = F)


# Table A.17
# Female Interviewer Controls
mod_interv_t = glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + nidaa + ennahda + intrvwr_f, family = binomial(), data=tun)
mod_interv_m <- glm(quota_dummy ~ gov_index + elec_free + female + gender_index_soc + islamist + edu + urban + household_fin + married + age + pjd + pam + intrvwr_f, family = binomial(), data=mor)
stargazer(mod_interv_t, mod_interv_m,
          table.placement="H", omit.stat=c("f", "ser", "rsq", "adj.rsq"),
          dep.var.labels.include = F,
          covariate.labels = c("Government \nPerformance", "Trust in \nElections", "Female", "Gender \nEgalitarianism", "Islamist", "Education", "Urban", "Household \nFinances", "Married", "Age", "Nidaa \nTounes", "Ennahda", "PJD", "PAM", "Female \nInterviewer"), no.space = F)
